# Load packages
library(haven)
library(tidyverse)
library(nleqslv)
library(data.table)
library(BB)

require("setRNG")
old.seed <- setRNG(list(kind="Mersenne-Twister", normal.kind="Inversion",
                        seed=123))

source(paste0(CODE,"parameters.R"))

# No spillovers exact hat estimation
counter<-function(hat_KAPPA,hat_KAPPA_s){
  
  source(paste0(CODE,"parameters.R"))
  df<-read_csv(paste(DATA, "clean_data4.csv", sep=""))
  
  areas<- unique(df$sub_area_o)
  df_temp<- df %>% group_by(sub_area_o,sub_area_d) %>% summarise(commute= sum(commute),commute_w= sum(commute_w),leisure= sum(leisure),leisure_w= sum(leisure_w))
  df_temp <- df_temp %>% group_by(sub_area_o) %>% mutate(share_commute = commute/sum(commute), share_commute_w = commute_w/sum(commute_w),share_leisure = leisure/sum(leisure), share_leisure_w = leisure_w/sum(leisure_w)) %>% ungroup()
  
  share_leisure<-df_temp$share_leisure
  share_leisure_w<-df_temp$share_leisure_w
  share_commute<-df_temp$share_commute
  share_commute_w<-df_temp$share_commute_w
  
  
  df_o2<- df_temp %>% group_by(sub_area_o) %>% summarise(residents = sum(commute),residents_w = sum(commute_w)) %>% ungroup() %>% mutate(share_resident=residents/sum(residents),share_resident_w=residents_w/sum(residents_w))
  df_o<-df %>% select(sub_area_o,share_resident,share_resident_w,residential_quant,avg_wage,avg_wage_w,res_fe_L,res_fe_L_w,res_fe_S,res_fe_S_w) %>% distinct()
  df_d<-df %>% select(sub_area_d,commercial_quant,office_quant,commercial_rent) %>% distinct() %>% group_by(sub_area_d) %>% summarise(commercial_sum=sum(commercial_quant*commercial_rent),office_sum=sum(office_quant*commercial_rent),commercial_quant=sum(commercial_quant),office_quant=sum(office_quant))
  
  share_resident<-df_o2$share_resident
  share_resident_w<-df_o2$share_resident_w
  
  
  share_N_p<- ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_p)/( (1-beta_N)*beta_T*beta_T_p) )/(1+ ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_p)/( (1-beta_N)*beta_T*beta_T_p) ))
  share_T_p<-1-share_N_p
  share_N_m<-((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_m)/( (1-beta_N)*beta_T*beta_T_m) )/(1+ ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_m)/( (1-beta_N)*beta_T*beta_T_m) ))
  share_T_m<-1-share_N_m
  
  residential_quant<-df_o$residential_quant
  office_quant<-df_d$office_quant
  commercial_quant<-df_d$commercial_quant
  
  BIGW<-gamma((epilson_l-1)/epilson_l)*exp(-df_o$res_fe_L)^(1/epilson_l)
  BIGS<-gamma((epilson_s-1)/epilson_s)*exp(-df_o$res_fe_S)^(alpha_N_p/epilson_s)
  BIGW_w<-gamma((epilson_l_w-1)/epilson_l_w)*exp(-df_o$res_fe_L_w)^(1/epilson_l_w)
  BIGS_w<-gamma((epilson_s_w-1)/epilson_s_w)*exp(-df_o$res_fe_S_w)^(alpha_N_m/epilson_s_w)
  
  kappa_s<-kappa
  kappa_s_w<-kappa_w
  
  # indices
  num<-sqrt(nrow(df_temp))
  num_x<-nrow(df_temp)
  xlength<-(6*num)
  x<-rep(1,xlength)
  own_index<- (1:num)+num*(0:(num-1))
  
  df_d<-df %>% select(sub_area_d,share_resident,share_resident_w,share_commute,share_commute_w) %>% group_by(sub_area_d) %>% summarise(total_share=sum(share_resident*share_commute),total_share_w=sum(share_resident_w*share_commute_w))
  df_d$temp<- df_d$total_share*share_N_p - df_d$total_share_w*share_N_m
  
  broydt <- function(x) {
    print(x)
    f<-rep(NA, length(x))
    
    hat_Q<-x[1:num]
    hat_w_N<-x[(num+1):(2*num)] 
    hat_w_T<-x[(2*num+1):(3*num)] 
    hat_w_w_T<-x[(3*num+1):(4*num)]
    hat_w_w_N<-x[(4*num+1):(5*num)]
    hat_P_N<-x[(5*num+1):(6*num)]
    
    hat_B<-1
    hat_B_w<-1
    hat_A_N<-1
    hat_A_T<-1
    
    hat_q=hat_Q
    
    temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute=share_commute,df_hat_w_N=rep(hat_w_N,num),df_hat_w_T=rep(hat_w_T,num),df_share_N_p=rep(share_N_p,num),df_share_T_p=rep(share_T_p,num),df_hat_KAPPA=hat_KAPPA)
    temp<- temp %>% group_by(origin) %>% mutate(out_N=(df_hat_w_N^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)),out_T=(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)), out2=  sum( df_commute*df_share_T_p*(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) + df_commute*df_share_N_p*(df_hat_w_N^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) ) )
    hat_share_commute_N<- temp$out_N/temp$out2
    hat_share_commute_T<-  temp$out_T/temp$out2 
    
    temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute_w=share_commute_w,df_hat_w_w_N=rep(hat_w_w_N,num),df_hat_w_w_T=rep(hat_w_w_T,num),df_share_N_m=rep(share_N_m,num),df_share_T_m=rep(share_T_m,num),df_hat_KAPPA=hat_KAPPA)
    temp<- temp %>% group_by(origin) %>% mutate(out_N=(df_hat_w_w_N^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)),out_T=(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)), out2=  sum( df_commute_w*df_share_T_m*(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) + df_commute_w*df_share_N_m*(df_hat_w_w_N^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) ) )
    
    hat_share_commute_w_N<-  temp$out_N/temp$out2
    hat_share_commute_w_T<-  temp$out_T/ temp$out2
    
    temp<-data_frame(origin=ceiling((1:num_x)/num),df_leisure=share_leisure,df_hat_P_N=rep(hat_P_N,num),df_hat_KAPPA=hat_KAPPA_s)
    temp<- temp %>% group_by(origin) %>% mutate(out=(df_hat_P_N^epilson_s)*(df_hat_KAPPA^(-kappa_s*epilson_s)) ,out2= sum(df_leisure*(df_hat_P_N^epilson_s)*(df_hat_KAPPA^(-kappa_s*epilson_s))))
    hat_share_leisure<- temp$out/temp$out2 
    
    temp<-data_frame(origin=ceiling((1:num_x)/num),df_leisure_w=share_leisure_w,df_hat_P_N=rep(hat_P_N,num),df_hat_KAPPA=hat_KAPPA_s)
    temp<- temp %>% group_by(origin) %>% mutate(out=(hat_P_N^epilson_s_w)*(df_hat_KAPPA^(-kappa_s_w*epilson_s_w)),out2= sum(df_leisure_w*(df_hat_P_N^epilson_s_w)*(df_hat_KAPPA^(-kappa_s_w*epilson_s_w))))
    hat_share_leisure_w<- temp$out/ temp$out2
    
    hat_BIGS <- (hat_KAPPA_s[own_index]^(-kappa_s*epilson_s))^(1/epilson_s) / ( hat_share_leisure[own_index]*hat_P_N^(epilson_s))^(alpha_N_p/epilson_s)
    hat_BIGS_w<-  (hat_KAPPA_s[own_index]^(-kappa_s_w*epilson_s_w) )^(1/epilson_s_w) / (hat_share_leisure_w[own_index]*hat_P_N^(epilson_s_w))^(alpha_N_m/epilson_s_w)
    hat_BIGW <- (hat_w_N^(epilson_l) * hat_KAPPA[own_index]^(-kappa*epilson_l))^(1/epilson_l) / (hat_share_commute_N[own_index])^(1/epilson_l) 
    hat_BIGW_w <- (hat_w_w_N^(epilson_l_w) *hat_KAPPA[own_index]^(-kappa_w*epilson_l_w))^(1/epilson_l_w) / (hat_share_commute_w_N[own_index])^(1/epilson_l_w)
    
    hat_share_resident<-  ((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p)/sum(share_resident*((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p))
    hat_share_resident_w<- ((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m) / sum(share_resident_w*((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m))
    
    temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_hat_commute_N=hat_share_commute_N,df_share_N_p=rep(share_N_p,num),df_share_T_p=rep(share_T_p,num),df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out_N=sum( (df_commute*df_share_N_p*df_hat_commute_N)^(1-(1/epilson_l_w))  *df_share_resident*df_hat_share_resident),out_T=sum( (df_commute*df_share_T_p*df_hat_commute_T)^(1-(1/epilson_l_w))  *df_share_resident*df_hat_share_resident),  out2_N=sum((df_commute*df_share_N_p)^(1-(1/epilson_l_w))  *df_share_resident),  out2_T=sum((df_commute*df_share_T_p)^(1-(1/epilson_l_w))  *df_share_resident))
    hat_N_N_p <- temp$out_N/temp$out2_N
    hat_N_T_p <- temp$out_T/temp$out2_T
    
    
    temp<-data_frame(dest=rep(1:num,num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_hat_share_commute_w_N=hat_share_commute_w_N,df_share_N_m=rep(share_N_m,num),df_share_T_m=rep(share_T_m,num),df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out_N=sum( (df_commute_w*df_share_N_m*df_hat_share_commute_w_N)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w),out_T=sum( (df_commute_w*df_share_T_m*df_hat_share_commute_w_T)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w),  out2_N=sum( (df_commute_w*df_share_N_m)^(1-(1/epilson_l_w)) *df_share_resident_w),  out2_T=sum((df_commute_w*df_share_T_m)^(1-(1/epilson_l_w)) *df_share_resident_w))
    hat_N_N_m <- temp$out_N/temp$out2_N
    hat_N_T_m <- temp$out_T/temp$out2_T
    
    hat_residential_quant<- (((alpha_H_m* total_pop_w*hat_share_resident_w*hat_BIGW_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*hat_share_resident*hat_BIGW*share_resident*BIGW))/(((alpha_H_m*total_pop_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*share_resident*BIGW))) )/hat_Q
    
    hat_office_quant <- ( hat_N_T_p*hat_w_T  )/hat_q
    hat_commercial_quant <- ( hat_N_N_p*hat_w_N   )/hat_q
    
    hat_total_quant<-((hat_residential_quant*residential_quant + hat_office_quant*office_quant + hat_commercial_quant*commercial_quant)/(residential_quant + office_quant + commercial_quant))
    
    f[1:num] <- hat_Q - hat_total_quant^((1-mu_H)/mu_H)

    f[(num+1):(2*num)] <- hat_w_N - (hat_A_N*hat_P_N/ (hat_q^(1-beta_N) * hat_w_w_N^(beta_N*beta_N_m) ) )^(1/(beta_N*beta_N_p))
    f[(2*num+1):(3*num)] <- hat_w_T - (hat_A_T/ (hat_q^(1-beta_T) * hat_w_w_T^(beta_T*beta_T_m) ) )^(1/(beta_T*beta_T_p))
    
    f[(3*num+1):(4*num)] <- hat_N_T_m*hat_w_w_T - ( hat_N_T_p*hat_w_T  )
    f[(4*num+1):(5*num)] <- hat_N_N_m*hat_w_w_N - ( hat_N_N_p*hat_w_N  )
    
    temp<-data_frame(dest=rep(1:num,num),df_leisure=share_leisure,df_hat_leisure=hat_share_leisure,df_leisure_w=share_leisure_w,df_hat_leisure_w=hat_share_leisure_w,df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident=rep(share_resident,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident=rep(hat_share_resident,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out=sum(alpha_N_p * total_pop *df_leisure*df_hat_leisure*df_share_resident*df_hat_share_resident*df_BIGW*df_hat_BIGW + alpha_N_m * total_pop_w*df_leisure_w*df_hat_leisure_w*df_share_resident_w*df_hat_share_resident_w*df_BIGW_w*df_hat_BIGW_w), out2=sum(alpha_N_p * total_pop*df_leisure*df_share_resident*df_BIGW + alpha_N_m * total_pop_w*df_leisure_w*df_share_resident_w*df_BIGW_w))

    f[(5*num+1):(6*num)] <- hat_P_N*(hat_A_N*( (hat_N_T_m^beta_N_m * hat_N_T_p^beta_N_p)^beta_N)* (hat_commercial_quant^(1-beta_N))) - ((temp$out)/(temp$out2) )
    
    return(f)
  }
  
  p0 <- rep(1,xlength)
  temp<-broydt(p0)
  out1<-nleqslv(p0,broydt,control=list(maxit=1000))
  x<-out1$x
  
  hat_Q<-x[1:num]
  hat_w_N<-x[(num+1):(2*num)] 
  hat_w_T<-x[(2*num+1):(3*num)] 
  hat_w_w_T<-x[(3*num+1):(4*num)]
  hat_w_w_N<-x[(4*num+1):(5*num)]
  hat_P_N<-x[(5*num+1):(6*num)]
  
  hat_B<-1
  hat_B_w<-1
  
  hat_A_N<-1
  hat_A_T<-1
  
  hat_q=hat_Q
  
  temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute=share_commute,df_hat_w_N=rep(hat_w_N,num),df_hat_w_T=rep(hat_w_T,num),df_share_N_p=rep(share_N_p,num),df_share_T_p=rep(share_T_p,num),df_hat_KAPPA=hat_KAPPA)
  temp<- temp %>% group_by(origin) %>% mutate(out_N=(df_hat_w_N^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)),out_T=(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)), out2=  sum( df_commute*df_share_T_p*(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) + df_commute*df_share_N_p*(df_hat_w_N^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) ) )
  hat_share_commute_N<- temp$out_N/temp$out2
  hat_share_commute_T<-  temp$out_T/temp$out2 
  
  temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute_w=share_commute_w,df_hat_w_w_N=rep(hat_w_w_N,num),df_hat_w_w_T=rep(hat_w_w_T,num),df_share_N_m=rep(share_N_m,num),df_share_T_m=rep(share_T_m,num),df_hat_KAPPA=hat_KAPPA)
  temp<- temp %>% group_by(origin) %>% mutate(out_N=(df_hat_w_w_N^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)),out_T=(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)), out2=  sum( df_commute_w*df_share_T_m*(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) + df_commute_w*df_share_N_m*(df_hat_w_w_N^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) ) )
  
  hat_share_commute_w_N<-  temp$out_N/temp$out2
  hat_share_commute_w_T<-  temp$out_T/ temp$out2
  
  temp<-data_frame(origin=ceiling((1:num_x)/num),df_leisure=share_leisure,df_hat_P_N=rep(hat_P_N,num),df_hat_KAPPA=hat_KAPPA_s)
  temp<- temp %>% group_by(origin) %>% mutate(out=(df_hat_P_N^epilson_s)*(df_hat_KAPPA^(-kappa_s*epilson_s)) ,out2= sum(df_leisure*(df_hat_P_N^epilson_s)*(df_hat_KAPPA^(-kappa_s*epilson_s))))
  hat_share_leisure<- temp$out/temp$out2 
  
  temp<-data_frame(origin=ceiling((1:num_x)/num),df_leisure_w=share_leisure_w,df_hat_P_N=rep(hat_P_N,num),df_hat_KAPPA=hat_KAPPA_s)
  temp<- temp %>% group_by(origin) %>% mutate(out=(hat_P_N^epilson_s_w)*(df_hat_KAPPA^(-kappa_s_w*epilson_s_w)),out2= sum(df_leisure_w*(df_hat_P_N^epilson_s_w)*(df_hat_KAPPA^(-kappa_s_w*epilson_s_w))))
  hat_share_leisure_w<- temp$out/ temp$out2
  
  hat_BIGS <- (hat_KAPPA_s[own_index]^(-kappa_s*epilson_s))^(1/epilson_s) / ( hat_share_leisure[own_index]*hat_P_N^(epilson_s))^(alpha_N_p/epilson_s)
  hat_BIGS_w<-  (hat_KAPPA_s[own_index]^(-kappa_s_w*epilson_s_w) )^(1/epilson_s_w) / (hat_share_leisure_w[own_index]*hat_P_N^(epilson_s_w))^(alpha_N_m/epilson_s_w)
  hat_BIGW <- (hat_w_N^(epilson_l) * hat_KAPPA[own_index]^(-kappa*epilson_l))^(1/epilson_l) / (hat_share_commute_N[own_index])^(1/epilson_l) 
  hat_BIGW_w <- (hat_w_w_N^(epilson_l_w) *hat_KAPPA[own_index]^(-kappa_w*epilson_l_w))^(1/epilson_l_w) / (hat_share_commute_w_N[own_index])^(1/epilson_l_w)
  
  hat_share_resident<-  ((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p)/sum(share_resident*((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p))
  hat_share_resident_w<- ((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m) / sum(share_resident_w*((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m))
  
  temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_hat_commute_N=hat_share_commute_N,df_share_N_p=rep(share_N_p,num),df_share_T_p=rep(share_T_p,num),df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num))
  temp<- temp %>% group_by(dest) %>% summarise(out_N=sum( (df_commute*df_share_N_p*df_hat_commute_N)^(1-(1/epilson_l))  *df_share_resident*df_hat_share_resident),out_T=sum( (df_commute*df_share_T_p*df_hat_commute_T)^(1-(1/epilson_l))  *df_share_resident*df_hat_share_resident),  out2_N=sum((df_commute*df_share_N_p)^(1-(1/epilson_l))  *df_share_resident),  out2_T=sum((df_commute*df_share_T_p)^(1-(1/epilson_l))  *df_share_resident))
  hat_N_N_p <- temp$out_N/temp$out2_N
  hat_N_T_p <- temp$out_T/temp$out2_T
  
  temp<-data_frame(dest=rep(1:num,num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_hat_share_commute_w_N=hat_share_commute_w_N,df_share_N_m=rep(share_N_m,num),df_share_T_m=rep(share_T_m,num),df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
  temp<- temp %>% group_by(dest) %>% summarise(out_N=sum( (df_commute_w*df_share_N_m*df_hat_share_commute_w_N)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w),out_T=sum( (df_commute_w*df_share_T_m*df_hat_share_commute_w_T)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w),  out2_N=sum( (df_commute_w*df_share_N_m)^(1-(1/epilson_l_w)) *df_share_resident_w),  out2_T=sum((df_commute_w*df_share_T_m)^(1-(1/epilson_l_w)) *df_share_resident_w))
  hat_N_N_m <- temp$out_N/temp$out2_N
  hat_N_T_m <- temp$out_T/temp$out2_T
  
  hat_residential_quant<- (((alpha_H_m* total_pop_w*hat_share_resident_w*hat_BIGW_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*hat_share_resident*hat_BIGW*share_resident*BIGW))/(((alpha_H_m*total_pop_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*share_resident*BIGW))) )/hat_Q
  
  hat_office_quant <- ( hat_N_T_p*hat_w_T  )/hat_q
  hat_commercial_quant <- ( hat_N_N_p*hat_w_N   )/hat_q
  
  hat_total_quant<-((hat_residential_quant*residential_quant + hat_office_quant*office_quant + hat_commercial_quant*commercial_quant)/(residential_quant + office_quant + commercial_quant))
  
  S<-sum(share_resident*hat_share_resident* (hat_BIGS)*BIGS)/sum(share_resident*BIGS)
  S_w<-sum(share_resident_w*hat_share_resident_w* (hat_BIGS_w)*BIGS_w)/sum(share_resident_w*BIGS_w)
  
  W<-sum(share_resident*hat_share_resident* (hat_BIGW)*BIGW)/sum(share_resident*BIGW)
  W_w<-sum(share_resident_w*hat_share_resident_w* (hat_BIGW_w)*BIGW_w)/sum(share_resident_w*BIGW_w)
  
  U<-sum(share_resident* (hat_BIGW^epilson_R_p) * (hat_BIGS^epilson_R_p) * hat_Q^(-epilson_R_p*alpha_H_p))^(1/epilson_R_p)
  U_w<-sum(share_resident_w* (hat_BIGW_w^epilson_R_m) * (hat_BIGS_w^epilson_R_m) * hat_Q^(-epilson_R_m*alpha_H_m))^(1/epilson_R_m)
  source(paste0(CODE,"parameters.R"))
  
  return( list(U,U_w,hat_share_commute_N,hat_share_commute_T,hat_share_commute_w_N,hat_share_commute_w_T,hat_share_leisure,hat_share_leisure_w,hat_share_resident,hat_share_resident_w,hat_commercial_quant,hat_office_quant,hat_Q,S,S_w,W,W_w))
}

# No spillovers exact hat estimation - only commuting
counter_c<-function(hat_KAPPA,hat_KAPPA_s){
  
  df<-read_csv(paste(DATA, "clean_data4.csv", sep=""))
  
  areas<- unique(df$sub_area_o)
  df_temp<- df %>% group_by(sub_area_o,sub_area_d) %>% summarise(commute= sum(commute),commute_w= sum(commute_w),leisure= sum(leisure),leisure_w= sum(leisure_w))
  df_temp <- df_temp %>% group_by(sub_area_o) %>% mutate(share_commute = commute/sum(commute), share_commute_w = commute_w/sum(commute_w),share_leisure = leisure/sum(leisure), share_leisure_w = leisure_w/sum(leisure_w)) %>% ungroup()
  
  share_leisure<-df_temp$share_leisure
  share_leisure_w<-df_temp$share_leisure_w
  share_commute<-df_temp$share_commute
  share_commute_w<-df_temp$share_commute_w

  df_o2<- df_temp %>% group_by(sub_area_o) %>% summarise(residents = sum(commute),residents_w = sum(commute_w)) %>% ungroup() %>% mutate(share_resident=residents/sum(residents),share_resident_w=residents_w/sum(residents_w))
  df_o<-df %>% select(sub_area_o,share_resident,share_resident_w,residential_quant,avg_wage,avg_wage_w,res_fe_L,res_fe_L_w,res_fe_S,res_fe_S_w) %>% distinct()
  df_d<-df %>% select(sub_area_d,commercial_quant,office_quant,commercial_rent) %>% distinct() %>% group_by(sub_area_d) %>% summarise(commercial_sum=sum(commercial_quant*commercial_rent),office_sum=sum(office_quant*commercial_rent),commercial_quant=sum(commercial_quant),office_quant=sum(office_quant))
  
  share_resident<-df_o2$share_resident
  share_resident_w<-df_o2$share_resident_w
  
  
  share_N_p<- ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_p)/( (1-beta_N)*beta_T*beta_T_p) )/(1+ ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_p)/( (1-beta_N)*beta_T*beta_T_p) ))
  share_T_p<-1-share_N_p
  share_N_m<-((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_m)/( (1-beta_N)*beta_T*beta_T_m) )/(1+ ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_m)/( (1-beta_N)*beta_T*beta_T_m) ))
  share_T_m<-1-share_N_m
  
  residential_quant<-df_o$residential_quant
  office_quant<-df_d$office_quant
  commercial_quant<-df_d$commercial_quant
  
  BIGW<-gamma((epilson_l-1)/epilson_l)*exp(-df_o$res_fe_L)^(1/epilson_l)
  BIGS<-gamma((epilson_s-1)/epilson_s)*exp(-df_o$res_fe_S)^(alpha_N_p/epilson_s)
  BIGW_w<-gamma((epilson_l_w-1)/epilson_l_w)*exp(-df_o$res_fe_L_w)^(1/epilson_l_w)
  BIGS_w<-gamma((epilson_s_w-1)/epilson_s_w)*exp(-df_o$res_fe_S_w)^(alpha_N_m/epilson_s_w)
  
  kappa_s<-kappa
  kappa_s_w<-kappa_w
  
  alpha_N_m<-0
  alpha_N_p<-0
  alpha_T_m<-1-alpha_H_m
  alpha_T_p<-1-alpha_H_p
  share_leisure<-rep(0,nrow(df_temp))
  share_leisure_w<-rep(0,nrow(df_temp))
  BIGS<-rep(1,sqrt(nrow(df_temp)))
  BIGS_w<-rep(1,sqrt(nrow(df_temp)))
  share_N_p<- 0
  share_T_p<-1-share_N_p
  share_N_m<-0
  share_T_m<-1-share_N_m
  
  # indices
  num<-sqrt(nrow(df_temp))
  num_x<-nrow(df_temp)
  xlength<-(3*num)
  x<-rep(1,xlength)
  
  own_index<- (1:num)+num*(0:(num-1))
  
  df_d<-df %>% select(sub_area_d,share_resident,share_resident_w,share_commute,share_commute_w) %>% group_by(sub_area_d) %>% summarise(total_share=sum(share_resident*share_commute),total_share_w=sum(share_resident_w*share_commute_w))
  df_d$temp<- df_d$total_share*share_N_p - df_d$total_share_w*share_N_m
  
  broydt <- function(x) {
    print(x)
    f<-rep(NA, length(x))
    
    hat_Q<-x[1:num]
    hat_w_T<-x[(num+1):(2*num)] 
    hat_w_w_T<-x[(2*num+1):(3*num)] 
    hat_B<-1
    hat_B_w<-1
    hat_A_T<-1
    
    hat_q=hat_Q
    
    temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute=share_commute,df_hat_w_T=rep(hat_w_T,num),df_hat_KAPPA=hat_KAPPA)
    temp<- temp %>% group_by(origin) %>% mutate(out_T=(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)), out2=  sum( df_commute*(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) ) )
    hat_share_commute_T<-  temp$out_T/temp$out2 
    
    temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute_w=share_commute_w,df_hat_w_w_T=rep(hat_w_w_T,num),df_hat_KAPPA=hat_KAPPA)
    temp<- temp %>% group_by(origin) %>% mutate(out_T=(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)), out2=  sum( df_commute_w*(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) ) )
    
    hat_share_commute_w_T<-  temp$out_T/ temp$out2
    hat_share_leisure<- rep(1,num_x)
    hat_share_leisure_w<- rep(1,num_x)

    hat_BIGS <- rep(1,num)
    hat_BIGS_w<-  rep(1,num)
    hat_BIGW <- (hat_w_T^(epilson_l) * hat_KAPPA[own_index]^(-kappa*epilson_l))^(1/epilson_l) / (hat_share_commute_T[own_index])^(1/epilson_l) 
    hat_BIGW_w <- (hat_w_w_T^(epilson_l_w) *hat_KAPPA[own_index]^(-kappa_w*epilson_l_w))^(1/epilson_l_w) / (hat_share_commute_w_T[own_index])^(1/epilson_l_w)

    hat_share_resident<-  ((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p)/sum(share_resident*((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p))
    
    hat_share_resident_w<- ((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m) / sum(share_resident_w*((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m))

    temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out_T=sum( (df_commute*df_hat_commute_T)^(1-(1/epilson_l))  *df_share_resident*df_hat_share_resident), out2_T=sum((df_commute)^(1-(1/epilson_l))  *df_share_resident))
    hat_N_T_p <- temp$out_T/temp$out2_T
    
    temp<-data_frame(dest=rep(1:num,num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out_T=sum( (df_commute_w*df_hat_share_commute_w_T)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w), out2_T=sum((df_commute_w)^(1-(1/epilson_l_w)) *df_share_resident_w))
    hat_N_T_m <- temp$out_T/temp$out2_T

    hat_residential_quant<- (((alpha_H_m* total_pop_w*hat_share_resident_w*hat_BIGW_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*hat_share_resident*hat_BIGW*share_resident*BIGW))/(((alpha_H_m*total_pop_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*share_resident*BIGW))) )/hat_Q

    hat_office_quant <- ( hat_N_T_p*hat_w_T  )/hat_q
    
    hat_total_quant<-((hat_residential_quant*residential_quant + hat_office_quant*(office_quant + commercial_quant))/(residential_quant + office_quant + commercial_quant))

    f[1:num] <- hat_Q - hat_total_quant^((1-mu_H)/mu_H)

    f[(num+1):(2*num)] <- hat_w_T - (hat_A_T/ (hat_q^(1-beta_T) * hat_w_w_T^(beta_T*beta_T_m) ) )^(1/(beta_T*beta_T_p))

    f[(2*num+1):(3*num)] <- hat_N_T_m*hat_w_w_T - ( hat_N_T_p*hat_w_T  )

    return(f)
  }
  
  p0 <- rep(1,xlength)
  temp<-broydt(p0)
  out1<-nleqslv(p0,broydt,control=list(maxit=1000))
  x<-out1$x
  
  hat_Q<-x[1:num]
  hat_w_T<-x[(num+1):(2*num)] 
  hat_w_w_T<-x[(2*num+1):(3*num)] 
  hat_B<-1
  hat_B_w<-1
  hat_A_T<-1

  hat_q=hat_Q
  
  temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute=share_commute,df_hat_w_T=rep(hat_w_T,num),df_hat_KAPPA=hat_KAPPA)
  temp<- temp %>% group_by(origin) %>% mutate(out_T=(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)), out2=  sum( df_commute*(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) ) )
  hat_share_commute_T<-  temp$out_T/temp$out2 
  
  temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute_w=share_commute_w,df_hat_w_w_T=rep(hat_w_w_T,num),df_hat_KAPPA=hat_KAPPA)
  temp<- temp %>% group_by(origin) %>% mutate(out_T=(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)), out2=  sum( df_commute_w*(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) ) )
  
  hat_share_commute_w_T<-  temp$out_T/ temp$out2
  
  hat_share_leisure<- rep(1,num_x)
  hat_share_leisure_w<- rep(1,num_x)

  hat_BIGS <- rep(1,num)
  hat_BIGS_w<-  rep(1,num)
  hat_BIGW <- (hat_w_T^(epilson_l) * hat_KAPPA[own_index]^(-kappa*epilson_l))^(1/epilson_l) / (hat_share_commute_T[own_index])^(1/epilson_l) 
  hat_BIGW_w <- (hat_w_w_T^(epilson_l_w) *hat_KAPPA[own_index]^(-kappa_w*epilson_l_w))^(1/epilson_l_w) / (hat_share_commute_w_T[own_index])^(1/epilson_l_w)

  hat_share_resident<-  ((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p)/sum(share_resident*((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p))
  hat_share_resident_w<- ((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m) / sum(share_resident_w*((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m))
  temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num))
  temp<- temp %>% group_by(dest) %>% summarise(out_T=sum( (df_commute*df_hat_commute_T)^(1-(1/epilson_l))  *df_share_resident*df_hat_share_resident), out2_T=sum((df_commute)^(1-(1/epilson_l))  *df_share_resident))
  hat_N_T_p <- temp$out_T/temp$out2_T
  
  temp<-data_frame(dest=rep(1:num,num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
  temp<- temp %>% group_by(dest) %>% summarise(out_T=sum( (df_commute_w*df_hat_share_commute_w_T)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w), out2_T=sum((df_commute_w)^(1-(1/epilson_l_w)) *df_share_resident_w))
  hat_N_T_m <- temp$out_T/temp$out2_T

  hat_residential_quant<- (((alpha_H_m* total_pop_w*hat_share_resident_w*hat_BIGW_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*hat_share_resident*hat_BIGW*share_resident*BIGW))/(((alpha_H_m*total_pop_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*share_resident*BIGW))) )/hat_Q
  hat_office_quant <- ( hat_N_T_p*hat_w_T  )/hat_q
  hat_total_quant<-((hat_residential_quant*residential_quant + hat_office_quant*(office_quant + commercial_quant))/(residential_quant + office_quant + commercial_quant))
  
  S<-sum(share_resident*hat_share_resident* (hat_BIGS)*BIGS)/sum(share_resident*BIGS)
  S_w<-sum(share_resident_w*hat_share_resident_w* (hat_BIGS_w)*BIGS_w)/sum(share_resident_w*BIGS_w)
  
  W<-sum(share_resident*hat_share_resident* (hat_BIGW)*BIGW)/sum(share_resident*BIGW)
  W_w<-sum(share_resident_w*hat_share_resident_w* (hat_BIGW_w)*BIGW_w)/sum(share_resident_w*BIGW_w)
  
  U<-sum(share_resident* (hat_BIGW^epilson_R_p) * (hat_BIGS^epilson_R_p) * hat_Q^(-epilson_R_p*alpha_H_p))^(1/epilson_R_p)
  U_w<-sum(share_resident_w* (hat_BIGW_w^epilson_R_m) * (hat_BIGS_w^epilson_R_m) * hat_Q^(-epilson_R_m*alpha_H_m))^(1/epilson_R_m)
  
  source(paste0(CODE,"parameters.R"))
  
  return( list(U,U_w,1,hat_share_commute_T,1,hat_share_commute_w_T,1,1,hat_share_resident,hat_share_resident_w,1,hat_office_quant,hat_Q,S,S_w,W,W_w))
}

# With spillovers exact hat estimation
counter_agg<-function(hat_KAPPA,hat_KAPPA_s){
  
  source(paste0(CODE,"parameters.R"))
  df<-read_csv(paste(DATA, "clean_data4.csv", sep=""))
  
  areas<- unique(df$sub_area_o)
  df_temp<- df %>% group_by(sub_area_o,sub_area_d) %>% summarise(commute= sum(commute),commute_w= sum(commute_w),leisure= sum(leisure),leisure_w= sum(leisure_w))
  df_temp <- df_temp %>% group_by(sub_area_o) %>% mutate(share_commute = commute/sum(commute), share_commute_w = commute_w/sum(commute_w),share_leisure = leisure/sum(leisure), share_leisure_w = leisure_w/sum(leisure_w)) %>% ungroup()
  
  share_leisure<-df_temp$share_leisure
  share_leisure_w<-df_temp$share_leisure_w
  share_commute<-df_temp$share_commute
  share_commute_w<-df_temp$share_commute_w
  
  df_o2<- df_temp %>% group_by(sub_area_o) %>% summarise(residents = sum(commute),residents_w = sum(commute_w)) %>% ungroup() %>% mutate(share_resident=residents/sum(residents),share_resident_w=residents_w/sum(residents_w))
  
  df_o<-df %>% select(sub_area_o,share_resident,share_resident_w,residential_quant,avg_wage,avg_wage_w,res_fe_L,res_fe_L_w,res_fe_S,res_fe_S_w) %>% distinct()
  
  df_d<-df %>% select(sub_area_d,commercial_quant,office_quant,commercial_rent) %>% distinct() %>% group_by(sub_area_d) %>% summarise(commercial_sum=sum(commercial_quant*commercial_rent),office_sum=sum(office_quant*commercial_rent),commercial_quant=sum(commercial_quant),office_quant=sum(office_quant))
  
  share_resident<-df_o2$share_resident
  share_resident_w<-df_o2$share_resident_w
  
  share_N_p<- ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_p)/( (1-beta_N)*beta_T*beta_T_p) )/(1+ ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_p)/( (1-beta_N)*beta_T*beta_T_p) ))
  share_T_p<-1-share_N_p
  share_N_m<-((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_m)/( (1-beta_N)*beta_T*beta_T_m) )/(1+ ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_m)/( (1-beta_N)*beta_T*beta_T_m) ))
  share_T_m<-1-share_N_m
  
  residential_quant<-df_o$residential_quant
  office_quant<-df_d$office_quant
  commercial_quant<-df_d$commercial_quant
  
  BIGW<-gamma((epilson_l-1)/epilson_l)*exp(-df_o$res_fe_L)^(1/epilson_l)
  BIGS<-gamma((epilson_s-1)/epilson_s)*exp(-df_o$res_fe_S)^(alpha_N_p/epilson_s)
  BIGW_w<-gamma((epilson_l_w-1)/epilson_l_w)*exp(-df_o$res_fe_L_w)^(1/epilson_l_w)
  BIGS_w<-gamma((epilson_s_w-1)/epilson_s_w)*exp(-df_o$res_fe_S_w)^(alpha_N_m/epilson_s_w)
  
  kappa_s<-kappa
  kappa_s_w<-kappa_w
  
  # indices
  num<-sqrt(nrow(df_temp))
  num_x<-nrow(df_temp)
  xlength<-(8*num)
  x<-rep(1,xlength)
  
  own_index<- (1:num)+num*(0:(num-1))
  
  df_d<-df %>% select(sub_area_d,share_resident,share_resident_w,share_commute,share_commute_w) %>% group_by(sub_area_d) %>% summarise(total_share=sum(share_resident*share_commute),total_share_w=sum(share_resident_w*share_commute_w))
  df_d$temp<- df_d$total_share*share_N_p - df_d$total_share_w*share_N_m

  broydt <- function(x) {
    print(x)
    f<-rep(NA, length(x))
    
    hat_Q<-x[1:num]
    hat_w_N<-x[(num+1):(2*num)] 
    hat_w_T<-x[(2*num+1):(3*num)] 
    hat_w_w_T<-x[(3*num+1):(4*num)]
    hat_w_w_N<-x[(4*num+1):(5*num)]
    hat_P_N<-x[(5*num+1):(6*num)]
    hat_B<-x[(6*num+1):(7*num)]
    hat_B_w<-x[(7*num+1):(8*num)]
    
    hat_q=hat_Q

    temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute=share_commute,df_hat_w_N=rep(hat_w_N,num),df_hat_w_T=rep(hat_w_T,num),df_share_N_p=rep(share_N_p,num),df_share_T_p=rep(share_T_p,num),df_hat_KAPPA=hat_KAPPA)
    temp<- temp %>% group_by(origin) %>% mutate(out_N=(df_hat_w_N^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)),out_T=(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)), out2=  sum( df_commute*df_share_T_p*(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) + df_commute*df_share_N_p*(df_hat_w_N^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) ) )
    hat_share_commute_N<- temp$out_N/temp$out2
    hat_share_commute_T<-  temp$out_T/temp$out2 
    
    temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute_w=share_commute_w,df_hat_w_w_N=rep(hat_w_w_N,num),df_hat_w_w_T=rep(hat_w_w_T,num),df_share_N_m=rep(share_N_m,num),df_share_T_m=rep(share_T_m,num),df_hat_KAPPA=hat_KAPPA)
    temp<- temp %>% group_by(origin) %>% mutate(out_N=(df_hat_w_w_N^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)),out_T=(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)), out2=  sum( df_commute_w*df_share_T_m*(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) + df_commute_w*df_share_N_m*(df_hat_w_w_N^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) ) )
    
    hat_share_commute_w_N<-  temp$out_N/temp$out2
    hat_share_commute_w_T<-  temp$out_T/ temp$out2

    temp<-data_frame(origin=ceiling((1:num_x)/num),df_leisure=share_leisure,df_hat_P_N=rep(hat_P_N,num),df_hat_KAPPA=hat_KAPPA_s)
    temp<- temp %>% group_by(origin) %>% mutate(out=(df_hat_P_N^epilson_s)*(df_hat_KAPPA^(-kappa_s*epilson_s)) ,out2= sum(df_leisure*(df_hat_P_N^epilson_s)*(df_hat_KAPPA^(-kappa_s*epilson_s))))
    hat_share_leisure<- temp$out/temp$out2 
    
    temp<-data_frame(origin=ceiling((1:num_x)/num),df_leisure_w=share_leisure_w,df_hat_P_N=rep(hat_P_N,num),df_hat_KAPPA=hat_KAPPA_s)
    temp<- temp %>% group_by(origin) %>% mutate(out=(hat_P_N^epilson_s_w)*(df_hat_KAPPA^(-kappa_s_w*epilson_s_w)),out2= sum(df_leisure_w*(df_hat_P_N^epilson_s_w)*(df_hat_KAPPA^(-kappa_s_w*epilson_s_w))))
    hat_share_leisure_w<- temp$out/ temp$out2
    
    hat_BIGS <- (hat_KAPPA_s[own_index]^(-kappa_s*epilson_s))^(1/epilson_s) / ( hat_share_leisure[own_index]*hat_P_N^(epilson_s))^(alpha_N_p/epilson_s)
    hat_BIGS_w<-  (hat_KAPPA_s[own_index]^(-kappa_s_w*epilson_s_w) )^(1/epilson_s_w) / (hat_share_leisure_w[own_index]*hat_P_N^(epilson_s_w))^(alpha_N_m/epilson_s_w)
    hat_BIGW <- (hat_w_N^(epilson_l) * hat_KAPPA[own_index]^(-kappa*epilson_l))^(1/epilson_l) / (hat_share_commute_N[own_index])^(1/epilson_l) 
    hat_BIGW_w <- (hat_w_w_N^(epilson_l_w) *hat_KAPPA[own_index]^(-kappa_w*epilson_l_w))^(1/epilson_l_w) / (hat_share_commute_w_N[own_index])^(1/epilson_l_w)
    
    hat_share_resident<-  ((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p)/sum(share_resident*((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p))
    hat_share_resident_w<- ((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m) / sum(share_resident_w*((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m))
    
    temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_hat_commute_N=hat_share_commute_N,df_share_N_p=rep(share_N_p,num),df_share_T_p=rep(share_T_p,num),df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out_N=sum( (df_commute*df_share_N_p*df_hat_commute_N)^(1-(1/epilson_l_w))  *df_share_resident*df_hat_share_resident),out_T=sum( (df_commute*df_share_T_p*df_hat_commute_T)^(1-(1/epilson_l_w))  *df_share_resident*df_hat_share_resident),  out2_N=sum((df_commute*df_share_N_p)^(1-(1/epilson_l_w))  *df_share_resident),  out2_T=sum((df_commute*df_share_T_p)^(1-(1/epilson_l_w))  *df_share_resident))
    hat_N_N_p <- temp$out_N/temp$out2_N
    hat_N_T_p <- temp$out_T/temp$out2_T
    
    
    temp<-data_frame(dest=rep(1:num,num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_hat_share_commute_w_N=hat_share_commute_w_N,df_share_N_m=rep(share_N_m,num),df_share_T_m=rep(share_T_m,num),df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out_N=sum( (df_commute_w*df_share_N_m*df_hat_share_commute_w_N)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w),out_T=sum( (df_commute_w*df_share_T_m*df_hat_share_commute_w_T)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w),  out2_N=sum( (df_commute_w*df_share_N_m)^(1-(1/epilson_l_w)) *df_share_resident_w),  out2_T=sum((df_commute_w*df_share_T_m)^(1-(1/epilson_l_w)) *df_share_resident_w))
    hat_N_N_m <- temp$out_N/temp$out2_N
    hat_N_T_m <- temp$out_T/temp$out2_T
    
    temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_hat_commute_N=hat_share_commute_N,df_share_N_p=rep(share_N_p,num),df_share_T_p=rep(share_T_p,num),df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_hat_share_commute_w_N=hat_share_commute_w_N,df_share_N_m=rep(share_N_m,num),df_share_T_m=rep(share_T_m,num),df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out=sum( (df_commute_w*df_share_N_m*df_hat_share_commute_w_N *df_share_resident_w*df_hat_share_resident_w) + (df_commute_w*df_share_T_m*df_hat_share_commute_w_T *df_share_resident_w*df_hat_share_resident_w) + (df_commute*df_share_N_p*df_hat_commute_N *df_share_resident*df_hat_share_resident) + (df_commute*df_share_T_p*df_hat_commute_T *df_share_resident*df_hat_share_resident) ),  out2=sum((df_commute_w*df_share_N_m *df_share_resident_w) + (df_commute_w*df_share_T_m *df_share_resident_w) + (df_commute*df_share_N_p *df_share_resident) + (df_commute*df_share_T_p *df_share_resident) ))
    hat_A_N<- (temp$out/temp$out2)^mu_A
    hat_A_T<- (temp$out/temp$out2)^mu_A

    hat_residential_quant<- (((alpha_H_m* total_pop_w*hat_share_resident_w*hat_BIGW_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*hat_share_resident*hat_BIGW*share_resident*BIGW))/(((alpha_H_m*total_pop_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*share_resident*BIGW))) )/hat_Q

    hat_office_quant <- ( hat_N_T_p*hat_w_T  )/hat_q
    hat_commercial_quant <- ( hat_N_N_p*hat_w_N   )/hat_q
    
    hat_total_quant<-((hat_residential_quant*residential_quant + hat_office_quant*office_quant + hat_commercial_quant*commercial_quant)/(residential_quant + office_quant + commercial_quant))

    f[1:num] <- hat_Q - hat_total_quant^((1-mu_H)/mu_H)
    f[(num+1):(2*num)] <- hat_w_N - (hat_A_N*hat_P_N/ (hat_q^(1-beta_N) * hat_w_w_N^(beta_N*beta_N_m) ) )^(1/(beta_N*beta_N_p))
    f[(2*num+1):(3*num)] <- hat_w_T - (hat_A_T/ (hat_q^(1-beta_T) * hat_w_w_T^(beta_T*beta_T_m) ) )^(1/(beta_T*beta_T_p))
    f[(3*num+1):(4*num)] <- hat_N_T_m*hat_w_w_T - ( hat_N_T_p*hat_w_T  )
    f[(4*num+1):(5*num)] <- hat_N_N_m*hat_w_w_N - ( hat_N_N_p*hat_w_N  )

    temp<-data_frame(dest=rep(1:num,num),df_leisure=share_leisure,df_hat_leisure=hat_share_leisure,df_leisure_w=share_leisure_w,df_hat_leisure_w=hat_share_leisure_w,df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident=rep(share_resident,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident=rep(hat_share_resident,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out=sum(alpha_N_p * total_pop *df_leisure*df_hat_leisure*df_share_resident*df_hat_share_resident*df_BIGW*df_hat_BIGW + alpha_N_m * total_pop_w*df_leisure_w*df_hat_leisure_w*df_share_resident_w*df_hat_share_resident_w*df_BIGW_w*df_hat_BIGW_w), out2=sum(alpha_N_p * total_pop*df_leisure*df_share_resident*df_BIGW + alpha_N_m * total_pop_w*df_leisure_w*df_share_resident_w*df_BIGW_w))
    
    f[(5*num+1):(6*num)] <- hat_P_N*(hat_A_N*( (hat_N_T_m^beta_N_m * hat_N_T_p^beta_N_p)^beta_N)* (hat_commercial_quant^(1-beta_N))) - ((temp$out)/(temp$out2) )
    f[(6*num+1):(7*num)]<- hat_B - (hat_share_resident/hat_share_resident_w)^mu_U_p
    f[(7*num+1):(8*num)]<- hat_B_w - (hat_share_resident_w/hat_share_resident)^mu_U_m
    
    return(f)
  }
  
  p0 <- rep(1,xlength)
  temp<-broydt(p0)
  out1<-nleqslv(p0,broydt,control=list(maxit=1000))

  x<-out1$x
  
  hat_Q<-x[1:num]
  hat_w_N<-x[(num+1):(2*num)] 
  hat_w_T<-x[(2*num+1):(3*num)] 
  hat_w_w_T<-x[(3*num+1):(4*num)]
  hat_w_w_N<-x[(4*num+1):(5*num)]
  hat_P_N<-x[(5*num+1):(6*num)]
  hat_B<-x[(6*num+1):(7*num)]
  hat_B_w<-x[(7*num+1):(8*num)]
  
  hat_q=hat_Q

  temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute=share_commute,df_hat_w_N=rep(hat_w_N,num),df_hat_w_T=rep(hat_w_T,num),df_share_N_p=rep(share_N_p,num),df_share_T_p=rep(share_T_p,num),df_hat_KAPPA=hat_KAPPA)
  temp<- temp %>% group_by(origin) %>% mutate(out_N=(df_hat_w_N^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)),out_T=(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)), out2=  sum( df_commute*df_share_T_p*(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) + df_commute*df_share_N_p*(df_hat_w_N^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) ) )
  hat_share_commute_N<- temp$out_N/temp$out2
  hat_share_commute_T<-  temp$out_T/temp$out2 
  
  temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute_w=share_commute_w,df_hat_w_w_N=rep(hat_w_w_N,num),df_hat_w_w_T=rep(hat_w_w_T,num),df_share_N_m=rep(share_N_m,num),df_share_T_m=rep(share_T_m,num),df_hat_KAPPA=hat_KAPPA)
  temp<- temp %>% group_by(origin) %>% mutate(out_N=(df_hat_w_w_N^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)),out_T=(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)), out2=  sum( df_commute_w*df_share_T_m*(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) + df_commute_w*df_share_N_m*(df_hat_w_w_N^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) ) )
  
  hat_share_commute_w_N<-  temp$out_N/temp$out2
  hat_share_commute_w_T<-  temp$out_T/ temp$out2
  
  temp<-data_frame(origin=ceiling((1:num_x)/num),df_leisure=share_leisure,df_hat_P_N=rep(hat_P_N,num),df_hat_KAPPA=hat_KAPPA_s)
  temp<- temp %>% group_by(origin) %>% mutate(out=(df_hat_P_N^epilson_s)*(df_hat_KAPPA^(-kappa_s*epilson_s)) ,out2= sum(df_leisure*(df_hat_P_N^epilson_s)*(df_hat_KAPPA^(-kappa_s*epilson_s))))
  hat_share_leisure<- temp$out/temp$out2 
  
  temp<-data_frame(origin=ceiling((1:num_x)/num),df_leisure_w=share_leisure_w,df_hat_P_N=rep(hat_P_N,num),df_hat_KAPPA=hat_KAPPA_s)
  temp<- temp %>% group_by(origin) %>% mutate(out=(hat_P_N^epilson_s_w)*(df_hat_KAPPA^(-kappa_s_w*epilson_s_w)),out2= sum(df_leisure_w*(df_hat_P_N^epilson_s_w)*(df_hat_KAPPA^(-kappa_s_w*epilson_s_w))))
  hat_share_leisure_w<- temp$out/ temp$out2
  
  hat_BIGS <- (hat_KAPPA_s[own_index]^(-kappa_s*epilson_s))^(1/epilson_s) / ( hat_share_leisure[own_index]*hat_P_N^(epilson_s))^(alpha_N_p/epilson_s)
  hat_BIGS_w<-  (hat_KAPPA_s[own_index]^(-kappa_s_w*epilson_s_w) )^(1/epilson_s_w) / (hat_share_leisure_w[own_index]*hat_P_N^(epilson_s_w))^(alpha_N_m/epilson_s_w)
  hat_BIGW <- (hat_w_N^(epilson_l) * hat_KAPPA[own_index]^(-kappa*epilson_l))^(1/epilson_l) / (hat_share_commute_N[own_index])^(1/epilson_l) 
  hat_BIGW_w <- (hat_w_w_N^(epilson_l_w) *hat_KAPPA[own_index]^(-kappa_w*epilson_l_w))^(1/epilson_l_w) / (hat_share_commute_w_N[own_index])^(1/epilson_l_w)
  
  hat_share_resident<-  ((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p)/sum(share_resident*((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p))
  hat_share_resident_w<- ((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m) / sum(share_resident_w*((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m))

  temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_hat_commute_N=hat_share_commute_N,df_share_N_p=rep(share_N_p,num),df_share_T_p=rep(share_T_p,num),df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num))
  temp<- temp %>% group_by(dest) %>% summarise(out_N=sum( (df_commute*df_share_N_p*df_hat_commute_N)^(1-(1/epilson_l_w))  *df_share_resident*df_hat_share_resident),out_T=sum( (df_commute*df_share_T_p*df_hat_commute_T)^(1-(1/epilson_l_w))  *df_share_resident*df_hat_share_resident),  out2_N=sum((df_commute*df_share_N_p)^(1-(1/epilson_l_w))  *df_share_resident),  out2_T=sum((df_commute*df_share_T_p)^(1-(1/epilson_l_w))  *df_share_resident))
  hat_N_N_p <- temp$out_N/temp$out2_N
  hat_N_T_p <- temp$out_T/temp$out2_T
  
  temp<-data_frame(dest=rep(1:num,num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_hat_share_commute_w_N=hat_share_commute_w_N,df_share_N_m=rep(share_N_m,num),df_share_T_m=rep(share_T_m,num),df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
  temp<- temp %>% group_by(dest) %>% summarise(out_N=sum( (df_commute_w*df_share_N_m*df_hat_share_commute_w_N)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w),out_T=sum( (df_commute_w*df_share_T_m*df_hat_share_commute_w_T)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w),  out2_N=sum( (df_commute_w*df_share_N_m)^(1-(1/epilson_l_w)) *df_share_resident_w),  out2_T=sum((df_commute_w*df_share_T_m)^(1-(1/epilson_l_w)) *df_share_resident_w))
  hat_N_N_m <- temp$out_N/temp$out2_N
  hat_N_T_m <- temp$out_T/temp$out2_T
  
  temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_hat_commute_N=hat_share_commute_N,df_share_N_p=rep(share_N_p,num),df_share_T_p=rep(share_T_p,num),df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_hat_share_commute_w_N=hat_share_commute_w_N,df_share_N_m=rep(share_N_m,num),df_share_T_m=rep(share_T_m,num),df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
  temp<- temp %>% group_by(dest) %>% summarise(out=sum( (df_commute_w*df_share_N_m*df_hat_share_commute_w_N *df_share_resident_w*df_hat_share_resident_w) + (df_commute_w*df_share_T_m*df_hat_share_commute_w_T *df_share_resident_w*df_hat_share_resident_w) + (df_commute*df_share_N_p*df_hat_commute_N *df_share_resident*df_hat_share_resident) + (df_commute*df_share_T_p*df_hat_commute_T *df_share_resident*df_hat_share_resident) ),  out2=sum((df_commute_w*df_share_N_m *df_share_resident_w) + (df_commute_w*df_share_T_m *df_share_resident_w) + (df_commute*df_share_N_p *df_share_resident) + (df_commute*df_share_T_p *df_share_resident) ))
  hat_A_N<- (temp$out/temp$out2)^mu_A
  hat_A_T<- (temp$out/temp$out2)^mu_A

  hat_residential_quant<- (((alpha_H_m* total_pop_w*hat_share_resident_w*hat_BIGW_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*hat_share_resident*hat_BIGW*share_resident*BIGW))/(((alpha_H_m*total_pop_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*share_resident*BIGW))) )/hat_Q
  
  hat_office_quant <- ( hat_N_T_p*hat_w_T  )/hat_q
  hat_commercial_quant <- ( hat_N_N_p*hat_w_N   )/hat_q
  
  hat_total_quant<-((hat_residential_quant*residential_quant + hat_office_quant*office_quant + hat_commercial_quant*commercial_quant)/(residential_quant + office_quant + commercial_quant))
  
  S<-sum(share_resident*hat_share_resident* (hat_BIGS)*BIGS)/sum(share_resident*BIGS)
  S_w<-sum(share_resident_w*hat_share_resident_w* (hat_BIGS_w)*BIGS_w)/sum(share_resident_w*BIGS_w)
  
  W<-sum(share_resident*hat_share_resident* (hat_BIGW)*BIGW)/sum(share_resident*BIGW)
  W_w<-sum(share_resident_w*hat_share_resident_w* (hat_BIGW_w)*BIGW_w)/sum(share_resident_w*BIGW_w)
  
  U<-sum(share_resident* (hat_B^epilson_R_p)*(hat_BIGW^epilson_R_p) * (hat_BIGS^epilson_R_p) * hat_Q^(-epilson_R_p*alpha_H_p))^(1/epilson_R_p)
  U_w<-sum(share_resident_w* (hat_B_w^epilson_R_m)* (hat_BIGW_w^epilson_R_m) * (hat_BIGS_w^epilson_R_m) * hat_Q^(-epilson_R_m*alpha_H_m))^(1/epilson_R_m)
  source(paste0(CODE,"parameters.R"))
  return( list(U,U_w,hat_share_commute_N,hat_share_commute_T,hat_share_commute_w_N,hat_share_commute_w_T,hat_share_leisure,hat_share_leisure_w,hat_share_resident,hat_share_resident_w,hat_commercial_quant,hat_office_quant,hat_Q,S,S_w,W,W_w))
}

# With spillovers exact hat estimation - only commuting
counter_c_agg<-function(hat_KAPPA,hat_KAPPA_s){
  
  df<-read_csv(paste(DATA, "clean_data4.csv", sep=""))
  areas<- unique(df$sub_area_o)
  df_temp<- df %>% group_by(sub_area_o,sub_area_d) %>% summarise(commute= sum(commute),commute_w= sum(commute_w),leisure= sum(leisure),leisure_w= sum(leisure_w))
  df_temp <- df_temp %>% group_by(sub_area_o) %>% mutate(share_commute = commute/sum(commute), share_commute_w = commute_w/sum(commute_w),share_leisure = leisure/sum(leisure), share_leisure_w = leisure_w/sum(leisure_w)) %>% ungroup()
  
  share_leisure<-df_temp$share_leisure
  share_leisure_w<-df_temp$share_leisure_w
  share_commute<-df_temp$share_commute
  share_commute_w<-df_temp$share_commute_w
  
  df_o2<- df_temp %>% group_by(sub_area_o) %>% summarise(residents = sum(commute),residents_w = sum(commute_w)) %>% ungroup() %>% mutate(share_resident=residents/sum(residents),share_resident_w=residents_w/sum(residents_w))
  df_o<-df %>% select(sub_area_o,share_resident,share_resident_w,residential_quant,avg_wage,avg_wage_w,res_fe_L,res_fe_L_w,res_fe_S,res_fe_S_w) %>% distinct()
  df_d<-df %>% select(sub_area_d,commercial_quant,office_quant,commercial_rent) %>% distinct() %>% group_by(sub_area_d) %>% summarise(commercial_sum=sum(commercial_quant*commercial_rent),office_sum=sum(office_quant*commercial_rent),commercial_quant=sum(commercial_quant),office_quant=sum(office_quant))
  
  share_resident<-df_o2$share_resident
  share_resident_w<-df_o2$share_resident_w
  
  share_N_p<- ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_p)/( (1-beta_N)*beta_T*beta_T_p) )/(1+ ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_p)/( (1-beta_N)*beta_T*beta_T_p) ))
  share_T_p<-1-share_N_p
  share_N_m<-((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_m)/( (1-beta_N)*beta_T*beta_T_m) )/(1+ ((df_d$commercial_sum/df_d$office_sum)*( (1-beta_T)*beta_N*beta_N_m)/( (1-beta_N)*beta_T*beta_T_m) ))
  share_T_m<-1-share_N_m
  
  residential_quant<-df_o$residential_quant
  office_quant<-df_d$office_quant
  commercial_quant<-df_d$commercial_quant
  
  BIGW<-gamma((epilson_l-1)/epilson_l)*exp(-df_o$res_fe_L)^(1/epilson_l)
  BIGS<-gamma((epilson_s-1)/epilson_s)*exp(-df_o$res_fe_S)^(alpha_N_p/epilson_s)
  BIGW_w<-gamma((epilson_l_w-1)/epilson_l_w)*exp(-df_o$res_fe_L_w)^(1/epilson_l_w)
  BIGS_w<-gamma((epilson_s_w-1)/epilson_s_w)*exp(-df_o$res_fe_S_w)^(alpha_N_m/epilson_s_w)
  
  kappa_s<-kappa
  kappa_s_w<-kappa_w
  
  alpha_N_m<-0
  alpha_N_p<-0
  alpha_T_m<-1-alpha_H_m
  alpha_T_p<-1-alpha_H_p
  share_leisure<-rep(0,nrow(df_temp))
  share_leisure_w<-rep(0,nrow(df_temp))
  BIGS<-rep(1,sqrt(nrow(df_temp)))
  BIGS_w<-rep(1,sqrt(nrow(df_temp)))
  share_N_p<- 0
  share_T_p<-1-share_N_p
  share_N_m<-0
  share_T_m<-1-share_N_m
  
  # indices
  num<-sqrt(nrow(df_temp))
  num_x<-nrow(df_temp)
  xlength<-(5*num)
  x<-rep(1,xlength)
  own_index<- (1:num)+num*(0:(num-1))
  
  df_d<-df %>% select(sub_area_d,share_resident,share_resident_w,share_commute,share_commute_w) %>% group_by(sub_area_d) %>% summarise(total_share=sum(share_resident*share_commute),total_share_w=sum(share_resident_w*share_commute_w))
  df_d$temp<- df_d$total_share*share_N_p - df_d$total_share_w*share_N_m
  
  broydt <- function(x) {
    print(x)
    f<-rep(NA, length(x))
    
    hat_Q<-x[1:num]
    hat_w_T<-x[(num+1):(2*num)] 
    hat_w_w_T<-x[(2*num+1):(3*num)] 
    
    hat_B<-x[(3*num+1):(4*num)] 
    hat_B_w<-x[(4*num+1):(5*num)] 
    
    hat_q=hat_Q

    temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute=share_commute,df_hat_w_T=rep(hat_w_T,num),df_hat_KAPPA=hat_KAPPA)
    temp<- temp %>% group_by(origin) %>% mutate(out_T=(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)), out2=  sum( df_commute*(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) ) )
    hat_share_commute_T<-  temp$out_T/temp$out2 
    
    temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute_w=share_commute_w,df_hat_w_w_T=rep(hat_w_w_T,num),df_hat_KAPPA=hat_KAPPA)
    temp<- temp %>% group_by(origin) %>% mutate(out_T=(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)), out2=  sum( df_commute_w*(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) ) )
    
    hat_share_commute_w_T<-  temp$out_T/ temp$out2
    hat_share_leisure<- rep(1,num_x)
    hat_share_leisure_w<- rep(1,num_x)

    hat_BIGS <- rep(1,num)
    hat_BIGS_w<-  rep(1,num)
    hat_BIGW <- (hat_w_T^(epilson_l) * hat_KAPPA[own_index]^(-kappa*epilson_l))^(1/epilson_l) / (hat_share_commute_T[own_index])^(1/epilson_l) 
    hat_BIGW_w <- (hat_w_w_T^(epilson_l_w) *hat_KAPPA[own_index]^(-kappa_w*epilson_l_w))^(1/epilson_l_w) / (hat_share_commute_w_T[own_index])^(1/epilson_l_w)
    
    hat_share_resident<-  ((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p)/sum(share_resident*((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p))
    hat_share_resident_w<- ((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m) / sum(share_resident_w*((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m))

    temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out_T=sum( (df_commute*df_hat_commute_T)^(1-(1/epilson_l))  *df_share_resident*df_hat_share_resident), out2_T=sum((df_commute)^(1-(1/epilson_l))  *df_share_resident))
    hat_N_T_p <- temp$out_T/temp$out2_T

    temp<-data_frame(dest=rep(1:num,num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out_T=sum( (df_commute_w*df_hat_share_commute_w_T)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w), out2_T=sum((df_commute_w)^(1-(1/epilson_l_w)) *df_share_resident_w))
    hat_N_T_m <- temp$out_T/temp$out2_T
    
    temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_hat_commute_N=hat_share_commute_N,df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_hat_share_commute_w_N=hat_share_commute_w_N,df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
    temp<- temp %>% group_by(dest) %>% summarise(out=sum( (df_commute_w*df_hat_share_commute_w_T *df_share_resident_w*df_hat_share_resident_w) + (df_commute*df_hat_commute_T *df_share_resident*df_hat_share_resident) ),  out2=sum((df_commute_w *df_share_resident_w) + (df_commute *df_share_resident) ))
    hat_A_T<- (temp$out/temp$out2)^mu_A

    hat_residential_quant<- (((alpha_H_m* total_pop_w*hat_share_resident_w*hat_BIGW_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*hat_share_resident*hat_BIGW*share_resident*BIGW))/(((alpha_H_m*total_pop_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*share_resident*BIGW))) )/hat_Q

    hat_office_quant <- ( hat_N_T_p*hat_w_T  )/hat_q
    hat_total_quant<-((hat_residential_quant*residential_quant + hat_office_quant*(office_quant + commercial_quant))/(residential_quant + office_quant + commercial_quant))

    f[1:num] <- hat_Q - hat_total_quant^((1-mu_H)/mu_H)
    
    f[(num+1):(2*num)] <- hat_w_T - (hat_A_T/ (hat_q^(1-beta_T) * hat_w_w_T^(beta_T*beta_T_m) ) )^(1/(beta_T*beta_T_p))
    
    f[(2*num+1):(3*num)] <- hat_N_T_m*hat_w_w_T - ( hat_N_T_p*hat_w_T  )
    
    f[(3*num+1):(4*num)]<- hat_B - (hat_share_resident/hat_share_resident_w)^mu_U_p
    f[(4*num+1):(5*num)]<- hat_B_w - (hat_share_resident_w/hat_share_resident)^mu_U_m
    
    return(f)
  }
  
  p0 <- rep(1,xlength)
  temp<-broydt(p0)
  out1<-nleqslv(p0,broydt,control=list(maxit=1000))
 
  x<-out1$x
  
  hat_Q<-x[1:num]
  hat_w_T<-x[(num+1):(2*num)] 
  hat_w_w_T<-x[(2*num+1):(3*num)] 
  
  hat_B<-x[(3*num+1):(4*num)] 
  hat_B_w<-x[(4*num+1):(5*num)] 
  
  hat_q=hat_Q
 
  temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute=share_commute,df_hat_w_T=rep(hat_w_T,num),df_hat_KAPPA=hat_KAPPA)
  temp<- temp %>% group_by(origin) %>% mutate(out_T=(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)), out2=  sum( df_commute*(df_hat_w_T^epilson_l)*(df_hat_KAPPA^(-kappa*epilson_l)) ) )
  hat_share_commute_T<-  temp$out_T/temp$out2 
  
  temp<-data_frame(origin=ceiling((1:num_x)/num),df_commute_w=share_commute_w,df_hat_w_w_T=rep(hat_w_w_T,num),df_hat_KAPPA=hat_KAPPA)
  temp<- temp %>% group_by(origin) %>% mutate(out_T=(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)), out2=  sum( df_commute_w*(df_hat_w_w_T^epilson_l_w)*(df_hat_KAPPA^(-kappa_w*epilson_l_w)) ) )
  
  hat_share_commute_w_T<-  temp$out_T/ temp$out2
  
  hat_share_leisure<- rep(1,num_x)
  
  hat_share_leisure_w<- rep(1,num_x)

  hat_BIGS <- rep(1,num)
  hat_BIGS_w<-  rep(1,num)
  hat_BIGW <- (hat_w_T^(epilson_l) * hat_KAPPA[own_index]^(-kappa*epilson_l))^(1/epilson_l) / (hat_share_commute_T[own_index])^(1/epilson_l) 
  hat_BIGW_w <- (hat_w_w_T^(epilson_l_w) *hat_KAPPA[own_index]^(-kappa_w*epilson_l_w))^(1/epilson_l_w) / (hat_share_commute_w_T[own_index])^(1/epilson_l_w)

  hat_share_resident<-  ((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p)/sum(share_resident*((hat_B*hat_BIGS*hat_BIGW*(hat_Q^(-alpha_H_p)))^epilson_R_p))
  
  hat_share_resident_w<- ((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m) / sum(share_resident_w*((hat_B_w*hat_BIGS_w*hat_BIGW_w*(hat_Q^(-alpha_H_m)))^epilson_R_m))
  
  temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num))
  temp<- temp %>% group_by(dest) %>% summarise(out_T=sum( (df_commute*df_hat_commute_T)^(1-(1/epilson_l))  *df_share_resident*df_hat_share_resident), out2_T=sum((df_commute)^(1-(1/epilson_l))  *df_share_resident))
  hat_N_T_p <- temp$out_T/temp$out2_T
  
  
  temp<-data_frame(dest=rep(1:num,num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
  temp<- temp %>% group_by(dest) %>% summarise(out_T=sum( (df_commute_w*df_hat_share_commute_w_T)^(1-(1/epilson_l_w)) *df_share_resident_w*df_hat_share_resident_w), out2_T=sum((df_commute_w)^(1-(1/epilson_l_w)) *df_share_resident_w))
  hat_N_T_m <- temp$out_T/temp$out2_T
  
  temp<-data_frame(dest=rep(1:num,num),df_commute=share_commute,df_hat_commute_T=hat_share_commute_T,df_hat_commute_N=hat_share_commute_N,df_BIGW=rep(BIGW,each=num),df_hat_BIGW=rep(hat_BIGW,each=num),df_share_resident=rep(share_resident,each=num),df_hat_share_resident=rep(hat_share_resident,each=num),df_commute_w=share_commute_w,df_hat_share_commute_w_T=hat_share_commute_w_T,df_hat_share_commute_w_N=hat_share_commute_w_N,df_BIGW_w=rep(BIGW_w,each=num),df_hat_BIGW_w=rep(hat_BIGW_w,each=num),df_share_resident_w=rep(share_resident_w,each=num),df_hat_share_resident_w=rep(hat_share_resident_w,each=num))
  temp<- temp %>% group_by(dest) %>% summarise(out=sum( (df_commute_w*df_hat_share_commute_w_T *df_share_resident_w*df_hat_share_resident_w) + (df_commute*df_hat_commute_T *df_share_resident*df_hat_share_resident) ),  out2=sum((df_commute_w *df_share_resident_w) + (df_commute *df_share_resident) ))
  hat_A_T<- (temp$out/temp$out2)^mu_A

  hat_residential_quant<- (((alpha_H_m* total_pop_w*hat_share_resident_w*hat_BIGW_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*hat_share_resident*hat_BIGW*share_resident*BIGW))/(((alpha_H_m*total_pop_w*share_resident_w*BIGW_w) + (alpha_H_p*total_pop*share_resident*BIGW))) )/hat_Q

  hat_office_quant <- ( hat_N_T_p*hat_w_T  )/hat_q
  
  hat_total_quant<-((hat_residential_quant*residential_quant + hat_office_quant*(office_quant + commercial_quant))/(residential_quant + office_quant + commercial_quant))
  
  S<-sum(share_resident*hat_share_resident* (hat_BIGS)*BIGS)/sum(share_resident*BIGS)
  S_w<-sum(share_resident_w*hat_share_resident_w* (hat_BIGS_w)*BIGS_w)/sum(share_resident_w*BIGS_w)
  
  W<-sum(share_resident*hat_share_resident* (hat_BIGW)*BIGW)/sum(share_resident*BIGW)
  W_w<-sum(share_resident_w*hat_share_resident_w* (hat_BIGW_w)*BIGW_w)/sum(share_resident_w*BIGW_w)
  
  U<-sum(share_resident* (hat_B^epilson_R_p) * (hat_BIGW^epilson_R_p) * (hat_BIGS^epilson_R_p) * hat_Q^(-epilson_R_p*alpha_H_p))^(1/epilson_R_p)
  U_w<-sum(share_resident_w* (hat_B_w^epilson_R_m)* (hat_B_w* hat_BIGW_w^epilson_R_m) * (hat_BIGS_w^epilson_R_m) * hat_Q^(-epilson_R_m*alpha_H_m))^(1/epilson_R_m)
  
  source(paste0(CODE,"parameters.R"))
  
  return( list(U,U_w,1,hat_share_commute_T,1,hat_share_commute_w_T,1,1,hat_share_resident,hat_share_resident_w,1,hat_office_quant,hat_Q,S,S_w,W,W_w))
}